\defmodule{MultinormalPCAGen}

Extends \class{MultinormalGen} for a \emph{multivariate normal} distribution
\cite{tJOH72a}, generated via the method of principal components analysis
(PCA) of the covariance matrix. The covariance matrix $\boldSigma$ is 
decomposed (by the constructor) as $\boldSigma = \bV\boldLambda\bV^{\tr}$ where
 $\bV$ is an orthogonal matrix and $\boldLambda$ is the diagonal matrix made up
of the eigenvalues of $\boldSigma$. $\bV^{\tr}$ is the transpose 
matrix of $\bV$. The eigenvalues are ordered from the
largest ($\lambda_1$) to the smallest ($\lambda_d$). The random multinormal
vector $\bX$ is generated via
\[
  \bX = \bmu + \bA \bZ,
\]
where $\bA = \bV\sqrt{\boldLambda}$, and $\bZ$ is a $d$-dimensional vector of
independent standard normal random variates. The decomposition method
uses the \texttt{SingularValueDecomposition} class in \texttt{colt}.
% The constructor specifies how these normal variates are generated.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\bigskip\hrule

\begin{code}
\begin{hide}
/*
 * Class:        MultinormalPCAGen
 * Description:  multivariate normal random variable generator
 * Environment:  Java
 * Software:     SSJ 
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       
 * @since

 * SSJ is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License (GPL) as published by the
 * Free Software Foundation, either version 3 of the License, or
 * any later version.

 * SSJ is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.

 * A copy of the GNU General Public License is available at
   <a href="http://www.gnu.org/licenses">GPL licence site</a>.
 */
\end{hide}
package umontreal.iro.lecuyer.randvarmulti;

   import cern.colt.matrix.DoubleMatrix2D;
   import cern.colt.matrix.linalg.SingularValueDecomposition;
\begin{hide}
import umontreal.iro.lecuyer.probdist.NormalDist;
import umontreal.iro.lecuyer.randvar.NormalGen;
import umontreal.iro.lecuyer.rng.RandomStream;
import cern.colt.matrix.impl.DenseDoubleMatrix2D;
\end{hide}

public class MultinormalPCAGen extends MultinormalGen\begin{hide} {
   private double[] lambda;

   private static SingularValueDecomposition getSvd (DoubleMatrix2D sigma) {
      return (new SingularValueDecomposition (sigma));
   }

   private DoubleMatrix2D decompPCA (SingularValueDecomposition svd) {
      DoubleMatrix2D D = svd.getS ();
      // Calculer la racine carree des valeurs propres
      for (int i = 0; i < D.rows(); ++i) {
         lambda[i] = D.getQuick (i, i);
         D.setQuick (i, i, Math.sqrt (D.getQuick (i, i)));
      }
      DoubleMatrix2D P = svd.getV();
      return P.zMult (D, null);
   }

   private void initL() {
      if (mu.length != sigma.rows() || mu.length != sigma.columns())
         throw new IllegalArgumentException
            ("Incompatible mean vector and covariance matrix");
      lambda = new double[mu.length];
      sqrtSigma = decompPCA (getSvd(sigma));
   }
\end{hide}\end{code}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection* {Constructors}
\begin{code}

   public MultinormalPCAGen (NormalGen gen1, double[] mu, double[][] sigma)\begin{hide} {
      super(gen1, mu, sigma);
      initL();
   }\end{hide}
\end{code}
\begin{tabb}   Equivalent to
 \method{MultinormalPCAGen}{(NormalGen, double[], DoubleMatrix2D)}\texttt{(gen1, mu, new DenseDoubleMatrix2D(sigma))}.
\end{tabb}
\begin{code}

   public MultinormalPCAGen (NormalGen gen1, double[] mu,
                             DoubleMatrix2D sigma)\begin{hide} {
      super(gen1, mu, sigma);
      initL();
   }\end{hide}
\end{code}
\begin{tabb} Constructs a multinormal generator with mean vector \texttt{mu}
 and covariance matrix \texttt{sigma}. The mean vector must have the same
 length as the dimensions of the covariance matrix, which must be symmetric
 and positive semi-definite.
 If any of the above conditions is violated, an exception is thrown.
 The vector $\bZ$ is generated by calling $d$ times the generator \texttt{gen1},
 which must be a \emph{standard normal} 1-dimensional generator.
\end{tabb}
\begin{htmlonly}
   \param{gen1}{the one-dimensional generator}
   \param{mu}{the mean vector.}
   \param{sigma}{the covariance matrix.}
   \exception{NullPointerException}{if any argument is \texttt{null}.}
   \exception{IllegalArgumentException}{if the length of the mean
    vector is incompatible with the dimensions of the covariance matrix.}
\end{htmlonly}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
\subsubsection* {Methods}
\begin{code}
   public static DoubleMatrix2D decompPCA (double[][] sigma) \begin{hide}  {
      return decompPCA (new DenseDoubleMatrix2D (sigma));
   }\end{hide}
\end{code}
\begin{tabb} Computes the decomposition \texttt{sigma} =
$\boldSigma = \bV\boldLambda\bV^{\tr}$. Returns $\bA = \bV\sqrt{\boldLambda}$.
\end{tabb}
\begin{code}

   public static DoubleMatrix2D decompPCA (DoubleMatrix2D sigma) \begin{hide}  {
      // L'objet SingularValueDecomposition permet de recuperer la matrice
      // des valeurs propres en ordre decroissant et celle des vecteurs propres de
      // sigma (pour une matrice symetrique et definie-positive seulement).

      SingularValueDecomposition sv = getSvd (sigma);
      // D contient les valeurs propres sur la diagonale
      DoubleMatrix2D D = sv.getS ();
      // Calculer la racine carree des valeurs propres
      for (int i = 0; i < D.rows(); ++i)
         D.setQuick (i, i, Math.sqrt (D.getQuick (i, i)));
      DoubleMatrix2D P = sv.getV();
      // Multiplier par la matrice orthogonale (ici P)
      return P.zMult (D, null);
   }\end{hide}
\end{code}
\begin{tabb} Computes the decomposition \texttt{sigma} = $\boldSigma =
 \bV\boldLambda\bV^{\tr}$. Returns $\bA = \bV\sqrt{\boldLambda}$.
\end{tabb}
\begin{code}

   public DoubleMatrix2D getPCADecompSigma()\begin{hide} {
      return sqrtSigma.copy();
   }\end{hide}
\end{code}
\begin{tabb} Returns the matrix $\bA = \bV\sqrt{\boldLambda}$ of this object.
\end{tabb}
\begin{htmlonly}
   \return{the PCA square root of the covariance matrix}
\end{htmlonly}
\begin{code}

   public static double[] getLambda (DoubleMatrix2D sigma)\begin{hide} {
      SingularValueDecomposition sv = getSvd (sigma);
      DoubleMatrix2D D = sv.getS ();
      double[] lambd = new double[D.rows()];
      for (int i = 0; i < D.rows(); ++i)
         lambd[i] = D.getQuick (i, i);
      return lambd;
   }\end{hide}
\end{code}
\begin{tabb} Computes and returns the eigenvalues of \texttt{sigma} in
 decreasing order.
\end{tabb}
\begin{code}

   public double[] getLambda()\begin{hide} {
      return lambda;
   }\end{hide}
\end{code}
\begin{tabb} Returns the eigenvalues of $\boldSigma$ in decreasing order.
\end{tabb}
\begin{code}

   public void setSigma (DoubleMatrix2D sigma)\begin{hide} {
      if (sigma.rows() != mu.length || sigma.columns() != mu.length)
         throw new IllegalArgumentException
            ("Invalid dimensions of covariance matrix");
      this.sigma.assign (sigma);
      this.sqrtSigma = decompPCA(getSvd (sigma));
   }\end{hide}
\end{code}
\begin{tabb} Sets the covariance matrix $\bSigma$ of this multinormal generator
 to \texttt{sigma} (and recomputes $\bA$).
\end{tabb}
\begin{htmlonly}
   \param{sigma}{the new covariance matrix.}
   \exception{IllegalArgumentException}{if \texttt{sigma} has
    incorrect dimensions.}
\end{htmlonly}
\begin{code}

   public static void nextPoint (NormalGen gen1, double[] mu,
                                 DoubleMatrix2D sigma, double[] p)\begin{hide} {
      if (gen1 == null)
         throw new NullPointerException ("gen1 is null");

      NormalDist dist = (NormalDist) gen1.getDistribution();
      if (dist.getMu() != 0.0)
         throw new IllegalArgumentException ("mu != 0");
      if (dist.getSigma() != 1.0)
         throw new IllegalArgumentException ("sigma != 1");

      if (mu.length != sigma.rows() ||
          mu.length != sigma.columns())
         throw new IllegalArgumentException
            ("Incompatible mean vector and covariance matrix dimensions");
      double[] temp = new double[mu.length];
      DoubleMatrix2D sqrtSigma = decompPCA(sigma);
      for (int i = 0; i < temp.length; i++) {
         temp[i] = gen1.nextDouble();
         if (temp[i] == Double.NEGATIVE_INFINITY)
            temp[i] = -MYINF;
         if (temp[i] == Double.POSITIVE_INFINITY)
            temp[i] = MYINF;
      }
      for (int i = 0; i < temp.length; i++) {
         p[i] = 0;
         for (int c = 0; c < temp.length; c++)
            p[i] += sqrtSigma.getQuick (i, c)*temp[c];
         p[i] += mu[i];
      }
   }\end{hide}
\end{code}
\begin{tabb} Generates a $d$-dimensional vector from the multinormal 
 distribution with mean vector \texttt{mu} and covariance matrix
 \texttt{sigma}, using the one-dimensional normal generator \texttt{gen1} to
 generate the coordinates of $\bZ$, and using the PCA decomposition of
 $\bSigma$. The resulting vector is put into \texttt{p}.
 Note that this static method will be very slow for large dimensions, because
 it recomputes the singular value decomposition at every call. It is therefore
 recommended to use a \texttt{MultinormalPCAGen} object instead,
 if the method is to be called more than once.
\end{tabb}
\begin{htmlonly}
   \param{p}{the array to be filled with the generated point.}
   \exception{IllegalArgumentException}{if the one-dimensional normal
    generator uses a normal distribution with $\mu$ not equal to 0, or
    $\sigma$ not equal to 1.}
   \exception{IllegalArgumentException}{if the length of the mean
    vector is different from the dimensions of the covariance matrix,
    or if the covariance matrix is not symmetric and positive-definite.}
   \exception{NullPointerException}{if any argument is \texttt{null}.}
\end{htmlonly}
\begin{code}

   public static void nextPoint (NormalGen gen1, double[] mu,
                                 double[][] sigma, double[] p)\begin{hide} {
      nextPoint(gen1, mu, new DenseDoubleMatrix2D(sigma), p);
   }\end{hide}
\end{code}
\begin{tabb}   Equivalent to
 \method{nextPoint}{(NormalGen, double[], DoubleMatrix2D, double[])}\texttt{(gen1, mu, new DenseDoubleMatrix2D(sigma), p)}.
\end{tabb}
\begin{code}

   public void nextPoint (double[] p)\begin{hide} {
      int n = mu.length;
      for (int i = 0; i < n; i++) {
         temp[i] = gen1.nextDouble();
         if (temp[i] == Double.NEGATIVE_INFINITY)
            temp[i] = -MYINF;
         if (temp[i] == Double.POSITIVE_INFINITY)
            temp[i] = MYINF;
      }
      for (int i = 0; i < n; i++) {
         p[i] = 0;
         for (int c = 0; c < n; c++)
            p[i] += sqrtSigma.getQuick (i, c)*temp[c];
         p[i] += mu[i];
      }
   }
}\end{hide}
\end{code}
\begin{tabb} Generates a point from this multinormal distribution. This is much
faster than the static method as it computes the singular value decomposition
matrix only once in the constructor.
\end{tabb}
\begin{htmlonly}
   \param{p}{the array to be filled with the generated point}
\end{htmlonly}
